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Abstract 

Quantization in the inversion layer and phase coherent transport are antici- 
pated to have significant impact on device performance in 'ballistic' nanoscale 
transistors. While the role of some quantum effects have been analyzed qual- 
itatively using simple one dimensional ballistic models, two dimensional (2D) 
quantum mechanical simulation is important for quantitative results. In this 
paper, we present a framework for 2D quantum mechanical simulation of 
a nanotransistor / Metal Oxide Field Effect Transistor (MOSFET). This 
framework consists of the non equilibrium Green's function equations solved 
self-consistently with Poisson's equation. Solution of this set of equations is 
computationally intensive. An efficient algorithm to calculate the quantum 
mechanical 2D electron density has been developed. The method presented 
is comprehensive in that treatment includes the three open boundary condi- 
tions, where the narrow channel region opens into physically broad source, 
drain and gate regions. Results are presented for (i) drain current versus 
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drain and gate voltages, (ii) comparison to results from Medici, and (iii) gate 
tunneling current, using 2D potential profiles. Methods to reduce the gate 
leakage current are also discussed based on simulation results. 
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Typeset using REVT^X 



I. INTRODUCTION 



MOSFETs with channel lengths in the tens of nanometer regime have recently been 
demonstrated by various research labs. EH! Design considerations to yield devices with desir- 
able characteristics have been explored in references |]-|8]. Device physics of these MOSFETs 
were analysed using simple quasi one dimensional modelsJHil The best modeling approach 
for design and analysis of nanoscale MOSFETs is presently unclear, though a straightfor- 
ward application of semiclassical methods that disregards quantum mechanical effects is 
generally accepted to be inadequate. Quantum mechanical modeling of MOSFETs with 
channel lengths in the tens of nanometers is important for many reasons: 

(i) MOSFETs with ultrathin oxide require an accurate treatment of current injection from 
source, drain and gate. Gate leakage is important because it places a lower limit on the 
OFF current. 

(ii) Ballistic flow of electrons across the channel becomes increasingly important as the 
channel length decreases. 

(hi) The location of the inversion layer changes from the source to the drain end, and its 
role in determining the C-V and I-V characterestics is most accurately included by a self- 
consistent solution of Poisson's equation and a quantum mechanical description to compute 
the charge density. 

(iv) Approximate theories of quantum effects included in semi-classical MOSFET modeling 
tools are desirable from practical considerations because semi-classical methods are numer- 
ically less expensive, and much of the empirical and semi-classical MOSFET physics devel- 
oped over the last few decades continues to hold true in many regions of a nanoscale MOS- 
FET. Examples of semiclassical methods that consider some quantum mechanical aspects 
are the density gradient and effective potential^ approaches, and quantum mechanical 
approximations used in the Medici package.0 Fully quantum mechanical simulations can 
play an important role in benchmarking such simulators. 

Central to quantum mechanical approaches to charge transport modeling is self- 
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consistent solution of a wave equation to describe the quantum mechanical transport, Pois- 
son's equation, and equations for statistics of the particle ensemble. In the absence of 
electron-electron and electron-phonon interactions (state of the scatterer does not change), 
the Landauer-Buttiker formalism0@ is applicable. In this formalism, the wave equation 
is Schrodinger's equation and the statistics is represented throughout the device by the 
Fermi-Dirac distribution of particles incident from the contacts (source, drain and gate). 
In the presence of electron-phonon interaction, the Wigner function (WF) and non equilib- 
rium Green's function (NEGF) formalisms are used. The NEGF approach has been quite 
successful in modeling steady state transport in a wide variety of one dimensional (ID) 
semiconductor structures.@@ 

A number of groups have started developing theory and simulation for fully quantum 
mechanical two dimensional simulation of MOSFETs using the: real space approach ,EHEI 
k-space approach,!! Wigner function approach,!! and non equilibrium Green's function ap- 
proachliiEU Others groups have taken a hybrid approach using the Monte Carlo method. 
The Monte Carlo approach, has been quite successful in describing scattering mechanisms 
in MOSFETs, in comparison to fully quantum mechanical approaches, and can also in- 
clude ballistic effects and the role of quantized energy levels in the MOSFET inversion layer 
in an approximate manner.&B Discussing the relative merits of various approaches and 
quantum-corrected drift-diffusion approaches is important. In fact, such a comparison of 
methods using standard device structures has been initiatedll but much work remains to be 
done in comparing and studying the suitability of different methods. Comparison of various 
methods is not the purpose of this paper. The purpose of this paper is to describe devel- 
opment of a particular approach, namely the NEGF approach, for numerical simulation of 
MOSFETs with two dimensional (2D) doping profiles and charge injection from the source, 
drain and gate contacts. 2D simulation significantly increases computational effort over the 
ID case. Non-uniform spatial grids are essential to limit the total number of grid points 
while at the same time resolving physical features. A new algorithm for efficient computation 
of electron density without complete solution of the system of equations is presented. The 
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computer code developed is used to calculate the drain and gate tunneling current in ultra 
short channel MOSFETs. Results from our approach and Medici are compared. The paper 
is organized as follows: formalism (section |II|), role of polysilicon gate depletion [111 A] , slopes 
of Id versus the gate (V g ) and drain (Vd) voltages (sections [111 A| - 1111 C| ), role of anisotropic 
effective mass (section |illD| ), role of gate tunneling current as a function of gate oxide 
thickness and gate length in determining the OFF current (section |IIIE[) . It is emphasized 
that the calculations presented include a self-consistent treatment of two dimensional gate 
oxide tunneling. Prior treatments of gate oxide tunneling in semi-classical 2D simulators 
incorporated ID models. 



II. FORMALISM 

A. The governing equations 

We consider Nf, independent valleys for the electrons within the effective mass approxi- 
mation. The Hamiltonian of valley b is 

h 2 , d I 1 d\ d I 1 d\ d / 1 d\ 

where (m b x ,m b y ,m b z ) are the components of the effective mass in valley b. The equation of 
motion for the retarded {G r ) and less-than (G < ) Green's functions are@'0'@ 

[E - H^nW^r^E) - J drTr hub ,{r x ,T,E)G r ¥M {r,T 2 ,E) = S h>b2 6(n - r 2 ) (2) 

and 

[E - H^rJjG^in, f 2 , E)-J drZl^in, f, E)G< M (f, f 2 , E) = 

J drT 1 <^r l) T,E)G a blM {T,T 2 ,E), (3) 

where G a is the advanced Green's function. In the above equations, the coordinate spans 
only the device (see Fig. [TJ). The influence of the semi- infinite regions of the source (S), drain 
(D) and polysilicon gate (P), and scattering mechanisms (electron-phonon) are included via 
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the self-energy terms Ejj lj6 / and We assume that charge is injected independently from 

the contact into each valley. Then, b2 c = T, bi c 5 bltb2 , where C represents the self-energy 
due to contacts. Finally, the hole bands are treated using the drift-diffusion model, which 
is expected to be a good approximation for n-MOSFETs. 

The electrostatic potential varies in the (x, y) plane, and the system is translationally 
invariant along the z-axis. So, all quantities A(f 1,7%, E) depend only on the difference 
coordinate z\ — z 2 . Using the relation 

A(n,f 2 ,E) = J d ^e ik ^-^A(xi,yi,x 2 ,y2,k z ,E) , (4) 

the equations of motion for G r and G < simplify to 

[E - ^ _ Hb {ri)\G r b {fi,r2, k z , E) - J dr E r b (fi, f, k z , E)G r b {r, f 2 , k z , E) = 5(n - r 2 ) (5) 

and 

[E ~ ~ HbftWbft'?* fc » E ) ~ J df ^ r 2 , k z , E) = 

J drE<(f 1 ,f,k z ,E)G a b (f,f 2 ,k z ,E), (6) 

where Z b = Z b b , and for the remainder of the paper, r — > {x,y). 

The density of states [iV(r, A; z , E 1 )] and charge density [p(r, E 1 )] are the sum of the 
contributions from the individual valleys: 

N(f, k z , E)=J2 Mr, k z , E) = --Im[G£(r, r, k z , E)\ (7) 
b 71 

p{r,k z ,E) = J2pb(r,k z ,E) = -iG<{r,r,k z ,E) . (8) 

b 

B. G r and Discretized matrix equations 

Self-consistent solution of the Green's function and Poisson's equations requires repeated 
computation of the non-equilibrium charge density. This computation is often the most time 
consuming part in modeling the electronic properties of devices. 
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The common procedure to evaluate the electron density uses the expression 

p b (r, k z , E) = -iGf(f, r, k z , E) 

= -iJ dndf 2 G r b (f, ri, k z , E)E<(n, r 2 , k z , E)G a b {f 2) r, k z , E), (9) 

where G r (fi, f 2 , k z , E) must be computed between all N x N y grid points and those grid points 
involving a non zero S a . The operation count required to solve for all elements of G r scales 
as (N x N y ) 3 , and so the use of Eq. (|) is expensive. We have developed a new recursive 
algorithm to compute the electron density in systems when the discretized version of the 
LHS of Eqs. (|5|) and (||) is block tridiagonal. This algorithm requires only the evaluation 
of the diagonal blocks of G r . The operation count of this algorithm scales as N x N y (or 
N x Ny) when the diagonal blocks correspond to lattice points in the x (or y) direction. We 
summarize the recursive algorithm to calculate G r (section |11 (A ) as it sets the stage for the 



new algorithm to compute G < (section [[I D|) . We stress that Poisson's equation only requires 
the diagonal elements of G < (Eq. (HD). The algorithm we develop in section [II C| however 
computes the diagonal blocks of G < . While this is much better than using Eq. (|9]) directly 
as discussed above, new algorithms to solve for only the diagonal elements with operation 
counts smaller than N x N y (or N x N y ) are very desirable. 
In matrix form, Eqs. ([5]) and @ are written as 

A'G r = A and (10) 
and A'G< = £<G a . (11) 

The self-energies due to the S, D and P are non zero only along the lines y = ys = yi, 
V — Ud = VN y and x = xp respectively (see Fig. [I]). The A' matrix has a dimension of N x N y 
and is ordered such that all grid points located at a particular y-coordinate correspond to 
its diagonal blocks. The notation adopted is that A', ^(i, i') refers to the off-diagonal entry 
corresponding to grid points (xi,yjA and (x'^yjA. The non zero elements of the diagonal 
blocks of the A' matrix are given by 

A' jd (i, i)=E' - Vij - T jd (i + 1, i) - T jd {i - 1, i) - T j+hj (i, i) - T^ ld (i, i) 
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-S r s (x 4 , Xi)6 jtl - S r D (x i; Xi)5 jjNy - £$>(%•, _ sr ( x i> Vi\ ( 12 ) 
• \'.,!>' ± 1, i)=Tj,j(i ± 1, i) - £s(z i± i, Xi)<J,-,i - S^,(xj ± i, Xi)S jtNy 

-^(x^^yfx^yj) (13) 
^'(^ i ')= -S s( a: <> x i')$j,i ~ ^r>( x h Xi')8j,N y , for %' ^ i, i ± 1 , (14) 
where E' = E — h 2 k 2 /2m z and V^j = V{xi 1 yj). The off-diagonal blocks are 

4j-±ij(i, *) = T j±1J (i,i) - Ep(yj,yj±i)6i tl 

A' j>jl (i,i , ) = 0,foTj'^j, j±l. (15) 

The non zero elements of the T matrix are defined by 

h 2 2 1 

r 7 -,-(i±l,i) = — 1 r !6 

2m ±x x i+ i - |x i± i - Xi\ 



Lj+li !J = t : 

7 2m ± »y i+1 -y i _ 1 | % -±i- % -| 



(17) 



where m ±:r = —. and m ±v = —. . Non zero elements of Sp(w,,«'), where 

j' 7^ j are neglected to ensure that A' is block tridiagonal (the algorithm to calculate G r and 
G K relies on the block tridiagonal form of A'). The A appearing in Eq. ( |T0|) corresponds to 
the delta function in Eq. (H). A is a diagonal matrix whose elements are given by 

4 

= 7 w \ • (18) 

{x i+1 -Xi_i){y i+ i - 



C. Recursive algorithm to calculate G r 

Pre-multiplying Eq. ( |10D by A -1 , 

A G r = I , (19) 

where matrix A is a symmetric matrix for both uniform and non uniform rectangular grids 
(Note that A' is symmetric only for a uniform grid). The recursive algorithm to calculate 
the diagonal blocks of the full Green's function is discussed now, using Dyson's equation for 
G r , and the left-connected Green's function as in referencesSS 
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(i) Dyson's equation for G r : The solution to 



A z,z A zz , \ I G' zz G' zzl \ _ / I O 
A z',z A z',z> ) \ G T „, „ G T „, ) \ O I 



(20) 



is 



G r = qtO + QrO UG r ^ 

= G r0 + G r UG r0 , (22) 



where, 



The advanced Green's function (G a ) is by definition related to G r by 

G a = G r\ = G a0 + QoOjj^a ^ 

= G a0 + G a U ] G a0 . (25) 



Eq. (pip is called Dyson's equation.li2rE!l 

(ii) Left-connected retarded Green's function : The left-connected (superscript L) retarded 
(superscript r) Green's function g rLq is defined by the first q blocks of Eq. ( |19"1) (includes 
the open boundaries attached to the lattice points via the self-energy) by 

A V .q,l:q 9^ = Iq,q, where, I q = h:q,V.q ■ (26) 

grLq+i j g (-[ggngQi j n a manner identical to g rLq except that the left-connected system is 
comprised of the first q + 1 blocks of Eq. (|19|). In terms of Eq. (p0|), the equation governing 
grLq+i f ]2 0ws by setting Z = 1 : q and Z' = g + 1. Using Dyson's equation [Eq. (pT|)1, we 
obtain 

rLg+l _ ( A _i_ /I ^ r ^' 



flg+l.g+l — (^9+1,9+1 + ^.g+l^^/^-g.g+l J • (27) 

Note that the last element is equal to the fully connected Green's function G r N N , which 
is the solution to Eq. [Tj|. 
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(iii) Full Green's function in terms of the left-connected Green's function : Consider Eq. 
([20|) such that A z ,z = A 1:qA . q , A z >, z > = A q+1:N)q+1 . N and A ZjZ ' = A l:qs+l . N . Noting that the 
only nonzero element of Ai :qtq+ i : N is A m+ i and using Eq. we obtain 

G\ q = g q L q + Q r q L q \A q ,q+iG r q+l q+1 A q+liq ) g q L q q (28) 
= C q + C q A g>q+1 G q+1>q . (29) 

Both G r and G q+1 are used in the algorithm for electron density, and so storing both sets 
of matrices will be useful. 

In view of the above equations, the algorithm to compute the diagonal blocks G q q is 
given by the following steps: 

• g^ 1 = a-, 1 . 



For q — 1, 2, N — 1, compute Eq. (B7 



• For q — N — 1, N — 2, 1, compute Eq. (f29|). Store G q+lq if memory permits for use 
in the algorithm for electron density. 

D. Recursive algorithm to calculate density (G < ) 

The discretized form of Eq. (D) is 

A'G< = S<G a , (30) 
where the dimension of the matrices involved are N = N x N y . Pre-multiplying by A -1 , 

AG< = S<G a , (31) 



where in Eq. ( |3"TD is equal to A -1 times the S < that appears in Eqs. (Q) and 

Following subsection |11 C] , the algorithm to calculate the electron density (diagonal el- 
ements of G < ) is discussed in terms of a Dyson's equation for G < and the left-connected 



g< L : 
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(i) Dyson's equation for G < : The solution to 

lz,z A z z , \ / G< z G< z , \ _ / B< z E< z , \ / G| z GJ ^, \ 
z'.z A z',z' J V G z',z G z',z' / V E z'.z E z',z' / V G z',z G z',z' / 

can be written as 

G < = G rOjj G < + G ^ s<G a ; (33) 

where G r0 and U have been defined in Eqs. (^), and G K and G a are readily identifiable 
from Eq. ©. Using G a = G a0 + G a0 WG a , Eq. (||) can be written as 

G< = G <0 + G <0 U^G a + G r0 UG < (34) 
= G <0 + G r UG <0 + G<U ] G a0 , (35) 
where G <0 = G r0 S < G a0 . (36) 

(ii) Left-connected g < : g <Lq is the counter part of g rLq , and is defined by the first q blocks 
ofEq. fl3j: 

g< L « = E< g)1:? . (37) 

g<Lq+i j g defined j n a manner identical to ^ <L|? except that the left-connected system is 
comprised of the first q + 1 blocks of Eq. (|3l"l). In terms of Eq. (^), the equation governing 
g<Lq+i f ii ows ^ setting Z = 1 : q and Z' = q + 1. Using the Dyson's equations for G r and 
G K , gq+iq+i can be recursively obtained (derivation is presented in Appendix A) as 



<Lq+l _ rLq+1 sp< , < a^9+i , ri^q+i sp< ai^q+i , r^g+i v < a-Lcj-l-i 

~~ yg+l,g+l L^g+l.g+l ' °9+lJ <" 39+1,9+1 ^9+1,9 1/9,9+1 ' #9+1,9 ^q,q+l #9+1,9+1' 

which can be written in a more intuitive form as 

^9+1,9+1 + °9+l + ^g+l,g 9^ A?,g+1 + ^9+1,9 dq^ ^q,q+l flg+l^+l) (39) 



aLg+1 I rLg+1 yi< aLq+1 . rLg+1 v-i< aLg+1 



<Lg+l _ rLg+1 
i/9+1,9+1 _ i/g+1,9+1 



where (75^ = Ag+i )g #^ 9 Ag g+1 . Eq. (|39| ) has the physical meaning that gq+i^+i has contri- 
butions due to four injection functions: (i) an effective self-energy due to the left-connected 
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structure that ends at q, which is represented by cr q+ x, (h) the diagonal self-energy com- 
ponent at grid point q + 1 that enters Eq. (|31~D, and (iii) the two off-diagonal self-energy 
components involving grid points q and q + 1. 

(iii) Full less-than Green's function in terms of left-connected Green's function : Consider 
Eq. Q32p such that A z = A 1:q)1 . q , A' z = A q+ i :N<q+ i :N and A z ,z> = A Uqtq+ i :N . Noting that the 
only nonzero element of Ai :qtq+ i :N is Aj^+i and using Eq. (|4|), we obtain 

^q,q = 9q,q 9 + 9q,q 9 ^9,g+l ^9+1,9 + 9 qt q+lA q+l q G q q + g q q A q ^ q+ \G q+l q . (40) 

Using Eq. (^), G q+1 can be written in terms of G< +1 +1 and other known Green's functions 
as 

n< _ „<0 \ n r A n <0 ±f r A n <L 1 J_ C 1< A^ n aL 1 (A~i\ 

^q+l,q ~ 9q+l,q ' ^ q+l^l^+^q+l^ ' ^9+1,9+1 J ' 1 ?+ 1 .9»«,g ' ^ q+l,q+l^q,q+l9 q,q ' V**-) 

Substituting Eq. (HJ) in Eq. (gD and using Eqs. (gl|) and (H), we obtain 



G^ q — 9^ q q + 9 q L q (^g.g+l^+l.q+l^g+l.g) + #<^ 9 Aj^+l^Vl-l^ + ^9,9+1^9+l,9#< 



<Lq 
'9,9 



+ 



$9^9+1^9+1,9^9,9 ^9,9^9,9+1^9+1,9 ' (^) 



where 



9q,q+l ~ 9q,q^q,q+l9q+l,q+l 

(43) 

9q+l,q = 9l+l,q+l^q+l,q9q q ■ (44) 



The terms inside the square brackets of Eq. (|42|) are Hermitian conjugates of each other. 

In view of the above equations, the algorithm to compute the diagonal blocks of G < is 
given by the following steps: 

. g< L1 = ggEfgg. 

• For q = 1,2, ...,N- 1, compute Eq. (0). 

• For q = N - 1, JV - 2, 1, compute Eqs. (E3) - 
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The current density flowing between two neighboring blocks q and q + 1 is given by 
J(q -> q + 1, k z , E)=J2 Jb{q 9 + 1, K E) 

b 

2e r i 

= ~r ^2 ^ r T qjq+ iG£: q q+1 (k z , E) - T q+ ^ q G^. q+l q {k z , E) , (45) 

where T has been defined in Eqs. ([TBI) and fll7| ). The current that has leaked into the gate 
between any two blocks p and q is 



Jgate = J2 J b(P ~> P + 1 , k z, E) - ^ J b (q -> + 1, fe„ £) 



(46) 



6 6 

and the total gate leakage current obtained by choosing p and g near the source and drain 
ends of the device. 

E. Expressions for Contact Self-energies (£5, and £p) 

Potential and doping profiles in the semi-infinite regions to the (a) left of 'S' and right of 
'D' are equal to the value at q = 1 and N y respectively (Fig. |l]). That is, they do not vary 
as a function of the y-coordinate, and (b) top of the 'P' is equal to the value of the top most 
grid line of 'P' (Fig. []]). That is, they are not a function of the x-coordinate. The retarded 
surface Green's functions of these semi-infinite regions are calculated from Eq. (|19D, when 
the matrices involved are semi-infinite. All diagonal sub-matrices of the A matrix are equal 
to A\ t \, Am Vj n v and Ap, and all first upper off-diagonal matrices of the A matrix are equal 
to A\ 2 , A Ny _i jNy and Ap_ ltP , in the source, drain and polysilicon regions respectively. We 
spell out the entire matrix for the source semi-infinite regions below: 



/ 


• 


• 



















/ . . . 




• 


• 


• 
















M . . . 







• 


• 


• 













M . . . 












Ai.i 


Ai.s 










M . . . 













A 2 ,l 


A x ,\ 


At, 2 







• • • 
















A 2 ,i 


Ai,i 


A X 


2 




V 

















A 2 ,i 


A X 


J 





9-3,-3 9-3,-2 9-3,-1 9-3,0 

9-2,-3 9-2,-2 9-2,-1 9-2,0 

9-1,-3 9-1,-2 9-1,-1 9-1,0 

90,-3 90,-2 90,-1 90,0 / 



\ /«oooooo\ 

• 

• 

/ 

7 

/0 

\00 00 / / 



(47) 



The surface Green's function of these regions can be obtained by using methods in matrix 
algebra that transform the two dimensional wire representing the semi-infinite contacts with 
N x grid points to N x one dimensional wires. 
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The self-energies due to the contacts are: 



Z r s (k z ,E) = Ai t0 g 0fi (k z , E)A Q ,i 

^d(^z,E) = A Ny<Ny+ ig Ny+ i <Ny+ i(k z , E)A Ny+ i jNy 

Z r P (k z ,E) = Apgp(k z ,E)Ap 

T,<{k z ,E) = -2iA lfi Im [g 0t0 (k z ,E)]A 0il f s (E) 

Z>n(k z ,E) = -2iA Ny:Ny+1 Im \g Ny+hNy+l (k z , E) 

£<(fc 2 , E) = -2iA P Im [g P (k z , E)) A P f P {E), 

where fi{E) is the Fermi factor in contact % G S, D, P. 



A 



Ny+l,NyfD(E) 



to 



[E xy - H b {fi)]G r b {f 1} f 2 , E xy ) - J drEKn, r, E xy )G r b (f, f 2 , E xy ) = 5(n 



r 2 



and 



[E xy - J ff fe (f 1 )]G < (f 1 , r 2 , Exy) - I drZKn, f, E xy )G<(r, f 2 , k z , E) 

J dr^{n,r,E xy )G a b {r,f 2 ,E X y) 



(48) 
(49) 
(50) 
(51) 
(52) 
(53) 



When S"(f*i, r-j, k z , E) depends only on E xy = E — -^f-, then Eqs. (|^) and (|^) simplify 



(54) 
(55) 



(56) 



While solving the equations, to keep the problem two dimensional, m z has to be independent 
of (x,y). So, we assume m z (Si0 2 ) = m z (Si). 



III. RESULTS AND DISCUSSION 

The steady state characteristics of MOSFETS that are of practical interest are the drive 
current, OFF current, slope of drain current versus drain voltage, and threshold voltage. 
In this section, we show that quantum mechanical simulations yield significantly different 
results from drift- diffusion based methods. These differences arise because of the following 
quantum mechanical features: 
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(i) polysilicon gate depletion in a manner opposite to the classical case, 

(ii) dependence of the resonant levels in the channel on the gate voltage, 

(iii) tunneling of charge across the gate oxide and from source to drain, 

(iv) quasi-ballistic flow of electrons. 

The MIT well-tempered 25 nm device structure!! is chosen for the purpose of discussion 
(MIT 25 nm device structure!! is hereafter referred to as MIT25). The method and computer 
code developed can however handle a wide variety of two dimensional structures with many 
terminals. We first compare the potential profiles from a constant mobility drift-diffusion 
solution and our quantum calculations at equilibrium. The motivation for this comparison 
results from the observation that the classical and quantum potential profiles should be in 
reasonable agreement, if the doping density is significantly higher than the electron and hole 
densities and the boundary conditions are the same. The doping profile of MIT25 meets 
this requirement in the channel region at small V g , and we verify that the potential profiles 
are in reasonable agreement at y = (see 'Ql flat band' and 'DD flat band' of Fig. [|). 
The legend 'flat band' refers to the potential at x = — t ox being fixed at the applied gate 
potential. 

An index of abbreviations used follows: 
Length Scales: t ox - oxide thickness, Lp - polysilicon gate thickness in x-direction, Lp - 
boundary of substrate region in x-direction, L y - Poisson's and NEGF equations are solved 
from —Ly/2 to +L y /2, L g - length of polysilicon gate region in y-direction. 
Models: Ql - quantum mechanical calculations using an isotropic effective mass, Q3 - 
quantum mechanical calculations using an anisotropic effective mass, DD - drift diffusion, 
Flat band - potential in the polysilicon gate region is held fixed from x = —(t ox + Lp) to 
x = —t ox at the bulk value, q-poly - potential in the gate polysilicon region is held fixed at 
x = —(t ox + Lp) at the bulk value, and the potential is computed quantum mechanically 
(self-consistently) for x > —(t ox + Lp). c-poly - classical treatment of gate polysilicon region, 
as in DD. 

Current and voltage: Id - drain current, I g - gate current, Vd - drain voltage, V g - gate 
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voltage. 

The values of constants assumed to obtain the numerical results of this section are as 
follows, unless otherwise noted: 

Electron effective mass of silicon: 0.3283 (isotropic), 0.19 and 0.98 (anisotropic), Electron 
effective mass of SiC^: m x = m y = 0.5 and same as silicon in m z direction, Hole effective 
mass of silicon is 0.49, band gap of silicon (SiC^): 1.12 eV (8.8 eV), energy barrier between 
the silicon and the oxide AEc=3.1 eV, dielectric constant of Si (SiC^) is e,si=11.9 (3.9) and 
kT = 0.02585 eV. 

A. Id versus V g - Effect of polysilicon depletion region 

The quantum mechanically calculated electron density near the SiC>2 barrier in the 
polysilicon region is smaller than the uniform background doping density. This is because 
the electron wavefunction is small close to the barrier. As a result, the conduction band in 
the polysilicon gate bends in a direction opposite to that computed semi-classically (compare 
x and triangle in Fig. ^).@Izl The band bending in the polysilicon gate plays a significant 
role in determining the threshold voltage and OFF current. To emphasize the importance 
of band bending, we plot the drain current versus gate voltage calculated with the gate 
polysilicon region treated as (i) 'flat band' and (ii) 'q-poly'. We find that the computed 
current is larger in (ii) because quantum mechanical depletion of electrons in the polysilicon 
gate region close to the oxide causes lowering of the potential in the channel. The Id versus 
V g curve shifts by approximately the an amount equal to the band bending in the polysili- 
con gate, in comparison to the flat band case. This band bending, which is measured from 
— (Lp + t ox ) to — t ox at equilibrium, is about 130 meV at the given doping density (Fig. ||). 
The influence of bandgap narrowing has been neglected here. It must be mentioned that 
the bandgap narrowing effect will tend to make the quantum mechanical contribution to the 
polysilicon band bending just discussed smaller. Future work to determine the influence of 
bandgap narrowing is necessary. 
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Computationally, a 2D treatment of the polysilicon gate region is expensive because of 
the additional grid points required. Note that matrix inversion depends on the cube of the 
matrix dimension. We point out that for highly doped polysilicon gate (in the absence of 
gate tunneling), a shift in the Id(V g ) curve from (i) by the equilibrium ID built-in potential 
does a reasonable job of reproducing the quantum mechanical result (see triangles in Fig. 
HD. This approximation becomes progressively poorer with increase in gate voltage, as can 
be seen from the figure. This is true espcecially for a smaller polysilicon doping density such 
as 1E20. 

B. Id versus V g - Comparison to Medici 

In the absence of gate tunneling and inelastic tunneling, the quantum mechanical current 

is 

I d = ^JdE T SD {E) [f s (E) - f D (E)] , (57) 

where Tsd is the transmission probability from source to drain, and fs and fjj are the 
Fermi-Dirac factors in the source and drain respectively. The total transmission (Fig. |J) 
is step-like with integer values at the plateaus in-spite of the complicated two dimensional 
electrostatics. In visual terms, the energies at which the steps turn on are determined by 
an effective 'subband dependent' source injection barrier, in contrast to the source injection 
barrier in drift-diffusion calculations.!^ This subband dependent source injection barrier is 
simply the maximum energy of the subband between source and drain due to quantization 
in the direction perpendicular to the gate plane (x-direction of Fig. 1). From a practical 
view point, the following two issues are important in ballistic MOSFETs: (a) typically, the 
total transmission assumes integer value at an energy slightly above the maximum in 2D 
density of states as shown in the inset of Fig. f|, and (b) the steps develop over 50 meV 
(twice the room temperature thermal energy). So, the shape of the steps is important in 
determining the value of current. Assuming a sharp step in total transmission with integer 
values in a calculation of current as in reference § is not quite accurate. 
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We compare the results from our quantum simulations with published results from 
quantum-corrected Medici.0 To compare the quantum and classical results, an estimate 
of the energy of the first subband minima (E r i) from Fig. |], and the location of the classi- 
cal barrier height (E^classical)) (Fig. |5[) are useful (Eb(classical) shown is obtained from 
constant mobility simulations using Prophet). The main features of this comparison are: 

(a) Subthreshold region: The slope d[log(Id)]/dV g is smaller in the quantum case when 
compared to Medici (Fig. ^|). Further, the current resulting from the simple intuitive 
expression 

/ = I q0 (58) 

matches the quantum result quite accurately. I q o is a prefactor chosen to reproduce the 
current at V g = in Fig. |6|. This match is rationalized by noting that for the values of gate 
biases considered, E r \ is well above the source Fermi energy and E r2 is many kT (thermal 
energy) above E r \. The difference in slope between the classical and quantum results can 
be understood from the slower variation of E r % in comparison to Eb(classical) as a function 
of V g (Fig. |H). We also find that the decrease of E r \ with increases in gate voltage is slower 
than the barrier height determined from the quantum potential profiles. This arises because 
(neglecting 2D effects) E rl is determined by a triangular well (whose apex is the conduction 
band) that becomes progressively narrower with increase in gate voltage. 

(b) Large gate biases: The drain current and slope d[log(Id)]/dV g are larger in the quantum 
case. The higher dId/dV g at large gate voltages in the quantum case can be understood 
from the fact that E r \ is above the Fermi level while Eb(classical) is below, at V g = IV (the 
quantum current is proportional to exp(— (E r \ — Ep)/kT)). The mobility model assumed 
in the classical case also plays a role in determining the slope. 

C. Id versus Vd 

The values of dld/dVd and drive current are important in MOSFET applications because 
they determine switching speeds, i Figure |7| compares the drain current versus drain voltage 
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for V g = and V g = IV. The drive current (V g = IV) calculated using Ql with the 
polysilicon region treated in the flat band and q-poly approximations is more than 100% 
and 200% larger than the results in reference |35|. dld/dVd in the linear region is up to three 
times larger in Ql. The subthreshold drain current is smaller in Ql. We however expect 
that with decreasing channel length, the sub threshold Id will become larger than the Medici 
results due to quantum mechanical tunnelingJll 

D. Isotropic versus anisotropic effective mass 

The primary influence of anisotropic effective mass is to influence the energy of the 
subbands in the inversion layer. Valleys with the largest effective mass perpendicular to 
the oxide (0.98m*) have subband energies that are smaller than the isotropic effective mass 
case. We see from the plot of transmission versus energy (Fig. |8]) that the valleys with 
(m x = 0.98m*, m y = m z = 0.19m*) have resonance levels that are more than 50meV lower 
in energy than the isotropic effective mass case. The corresponding subthreshold current 
(Fig. [5]) is a few hundred percent larger than the value obtained from the isotropic effective 
mass case. This follows by noting that the subthreshold current depends on exp(—E r \/kT). 
The drive current (Fig. |9|) from the anisotropic effective mass case is more than twenty five 
percent larger than the isotropic effective mass case. Note that for large gate voltages the 
dependence of current on E r \ is sub exponential. We are not aware of any calculations that 
compare the relative importance of the current carrying capacity of electrons in the three 
inequivalent valleys. We find that the valley with the largest m x (=0.98m*) carries 89.22 % 
and 79.77 % of the current at V g equal to and IV respectively (Vd = IV). Thus all three 
valleys are necessary for an accurate calculation of the ballistic current. 

E. Gate leakage current 

A major problem in MOSFETs with ultra thin oxides is that tunneling from gate to 
drain will determine the OFF current. The gate leakage current versus y is plotted for the 
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MIT25 device in Fig. [10]. Note that while we use a value of 3.0 for the dielectric constant of 
Si0 2 , a value of 3.9 does not change the qualitative conclusions. At V g = 0V and = IV, 
the main path for leakage current is from the polysilicon gate contact on top of the oxide 
to the highly doped (n + ) regions associated with the drain (Source Drain Extension, SDE) 
as shown in Fig. [T^ (a). At non zero V g , there is also an appreciable tunneling from the 
highly doped n + regions near the source to the polysilicon region on top of the gate (Fig. 
|10| (a)). For t ox = 1.5 nm, gate tunneling increases the OFF current by about two orders of 
magnitude, and for smaller oxide thicknesses, the gate leakage current is significantly larger. 

We propose that the gate leakage current can be reduced by a factor of 10-100 without 
significantly compromising the drive current. The drive current in these ultra small MOS- 
FETs is primarily determined by the source injection barrier, or more correctly as discussed 
earlier by the resonant level at the source injection barrier. So any changes that result in 
a reduction of the gate leakage current should not significantly alter the location of the 
resonant level at the source injection barrier (and hence the drive current). Two methods 
(without regard to fabrication issues) that help in this direction are discussed below: 

(i) Shorter or asymmetric polysilicon gate region: We propose that the gate leakage 
current can be significantly reduced by using shorter gate lengths. The main feature of the 
shorter gate lengths is a small overlap between the polysilicon gate and the n + region near 
the drain. This is pictorially represented in Figs. [TT] (a) and (b) with 'long' and 'short' 
gate lengths. To simulate the long and short gate lengths, we consider the doping profile of 



MIT25 with L g = 25 nm and 50 nm (gate length in reference |35[). The OFF current and 
gate leakage current are plotted in Fig. |12[ We see that the gate leakage current reduces 
by more than an order of magnitude, and the drive current is within two percent of the 
L g = 50 nm case, as desired (see inset of Fig. [12]). The spatial profile of gate leakage 
current for L g = 25 nm is shown in Fig. |H] (b). Though the gate leakage current reduces 
significantly, a drawback of this scheme is the requirement for very short (approximately 
equal to the distance between highly doped region near source and drain) polysilicon gate 
lengths. A polysilicon gate placed asymmetrically with respect to y=0 such that its overlap 
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with the n + regions near the drain is small, will also serve to reduce the OFF current without 
compromising the drive current. 

(ii) Graded oxide: The second proposal is to use a graded oxide, which is thinner close to 
the source end and thicker close to the drain end (Fig. [11] (c)). The thinner oxide near the 
source is not expected to alter the source injection barrier significantly, while the tunneling 
rate from gate to drain will be significantly smaller because of the thicker oxide in the drain- 
gate overlap region. We consider an oxide that is 1.5 nm thick for y < +10 nm and 2.5 nm 
for y > 11 nm, with the thickness varying linearly in between. The polysilicon gate lengths is 
50 nm. Comparison of this device to the original MIT25 with an uniform oxide and L g = 50 
nm show that while the gate leakage current decreases by one order of magnitude, the drive 
current decreases by only 30 %. Further optimization of this device structure could yield a 
larger drive current, while keeping the gate leakage current small. 

IV. CONCLUDING REMARKS 

A modeling framework and computer code to calculate properties of ballistic MOSFETs 
with open boundaries at the source, drain and gate contacts have been developed. This 
includes an algorithm to compute the electron density using the NEGF equations that avoids 
solving for the entire G r matrix even in the presence of non zero self energies throughout 
the device. Note that the simulations presented are 2D in nature and also involve self- 
consistency. As a result, they were numerically intensive and were typically performed on 
sixteen to sixty four processors of an SGI Origin machine. 

The main results of this study are: 

(a) Polysilicon gate depletion causes the conduction band close to the oxide interface to 
bend in a manner opposite to the semi-classical case (Fig. 2). This causes a substantial 
shift in the location of the conduction band bottom in the channel, which gives rise to drain 
currents that are different from the semiclassical case by one to two orders of magnitude. 
Performing quantum mechanical calculations with a flat polysilicon region, and then shifting 
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the gate voltage axis (in Id versus V g ) by the quantum mechanical built-in voltage shown in 
Fig. 2 results in an order of magnitude better agreement with results from a quantum me- 
chanical treatment of the polysilicon region. This built-in voltage can simply be determined 
by ID simulations or an analytical expression. In reality, treatment of discrete dopants in 
the polysilicon region will give rise to results that are in between the 'flat band' and 'q-poly' 
cases presented in the paper. 

A quantum mechanical treatment of the polysilicon gate region results in an OFF current 
(V g = V and Vd — 1 V) that is more than 35 times larger than the OFF current from a 
flat band treatment of polysilicon region and published results!! based on a sophisticated 
semiclassical simulator. 

(b) Resonant levels in the channel the from source to drain increase the effective source 
injection barrier for ballistic electrons. Further, even in the ballistic limit the transmission 
versus energy reaches integer values over an energy range that could be many times the 
thermal energy. Knowledge of the detailed shape of transmission versus energy is important 
to accurately determine the ballistic current. The precise shape of these transmission steps 
depends on the details of the channel to source and drain overlap regions and the resulting 2D 
potential profile. Assuming a sharp step-like increases in the total transmission is incorrect. 

The slope dld/dVd, whose importance was emphasized in reference |7| and the drive current 
(at V g = IV) are about 300% larger than reported in reference [35]. Further, inclusion of 
anisotropic effective mass in our calculation makes the quantum results deviate further from 
the semiclassical results as shown in Fig. |9|. 

(c) Tunneling of charge across the gate oxide can put a limit on the OFF current. Models 
of the tunnel current for thin oxide MOSFETs are important. We model the gate leakage 
current in two dimensions and show that significant reduction in the OFF current is possible 
without altering the drive current significantly. This is accomplished by changing either the 
gate length (Fig. [11] b) or by introducing a graded oxide (Figs. [11] (c)). 

(d) Quasi-ballistic flow of electrons causes the slope of d[log(Id)}/dV g to be larger than the 
values obtained from drift-diffusion methods using field dependent mobility models. 
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This paper dealt with the modeling of steady state properties of nanoscale MOSFETs in 
the ballistic limit. Future work in quantum mechanical simulation of MOSFETs that are of 
importance include: (i) treatment of scattering mechanisms such as interface roughness and 
electron-phonon scattering,il'0 (ii) treatment of discrete impurity dopants,0'0 (iii) switch- 
ing behavior of MOSFETs / time-dependent simulation@@ and (iv) noise characteristics of 
nano transistors. 
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VI. APPENDIX A 



Derivation of Eqs. (38) and (\3Q): 

Using Dyson's equation for G < [Eq. (B^)], we obtain 



<Lq+l _ rO a] n < L 1+ 1 i n r0 y< al,q+l rO y< al,q+L /^n 

i/q+l.g+l — i/g+l^+l^g+l.gi/ij.g+l ^ yq+l,q+l^q+l,q+iyq+l,q+l ^ yq+l,q+l^q+l,qyq,q+l ■ l dy J 

Using Eq. (0), g^q+i 1 can be expressed in terms of gq+it+i, the quantity we are solving 
for and known Green's functions as 

<Lq+l _ <0 , „<0 4t aLq+1 , n <Lq a\ aLq+1 rLq a <Lq+l / fif) N 

i/9,9+1 — i/9,9+1 + y^q+lAj+l.qy^q+l + i/ g>(? Aj^+iyq+l^+l + yg,g A g,9+li/g+l,g+l ■ l UU J 

Substituting Eq. (|60|) in Eq. (|59D , we obtain 

"r _ ^7 r0 A , n rL9 4 J « <L ' 3+1 

J yg+l.q+l^+l^yg.q ^9,9+lJ i/g+l.g+l ~~ 

„r0 y , < aLq+1 . rO y , < /1 ai 9+ 1 

y^+l^+l^g+l.g+iyg+l.g+l 1/9+1,9+1 9+l>9y9>9+l 

i„rO 4 [ <0 , <0 a\ aLq+1 , <Lq a\ aLq+1 1 / fi1 \ 

' y^+l^+l^+l,? Ly<j,cj+1 "i" yq^+l^q+l^yq^q+l {Jq,q ■ n q,q+iyq+l,q+l\ • l U± J 

Using Eq. (p7|) and g^q+i = y"™ 9 ^ 9 +iy"g+i ? +i> which follows from Eq. fl3~6|) , we obtain 



aLq+l . rO v- , < „aLq+l 



aLq+1 
9q+l,q+l 



<Lq+l _ rLq+1 IW , n^jt 

yg+1,9+1 — yg+l,?+l [^g+l.g+l ^ ^-q+hqilqa ^9,9+1 

, rLg+1 v < aLg+1 . rLg+l j „rLqsp< aLq+1 

-T-yq+l,q+l^q+l, q yq,q+l ' 9 q+l,q+l A q+hqy q,q ^ q,q+iy q+l,q+l • 



(62) 



Noting that y-^ij 1 = ^1^+1^9+1,9^^, Eq. (H) can be written as 



n <Lq+l _ rLq+1 

y q +i,q+i — y q +i, q +i 



^9+1,9+1 + a q+X 



aLq+1 , rLq+1 aLq+1 
i/9+1,9+1 1* 1/9+1,9+1^+1,^9,9+1 

1 / £ ?+ 1 V< „»^9+l 
"ryg+l,g ^q^+iyq+l^+li 



where, 

°g+l = Aq+l t qg^ q q A q q+l 



(63) 



(64) 
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Figure Captions: 

Fig. |1|: The equations are solved in a 2D non uniform spatial grid, with semi-infinite 
boundaries as shown. Each column q comrises the diagonal blocks of Eqs. ( |19D and (|3T|). 
The electrostatic potential is held fixed at the begining of the semi-infinte regions closest to 
the device. 

Fig. Potential profile at the y=0 slice of MIT25, calculated by four different methods. 
Note the qualitative difference of the 'Ql q-poly' case due to electron depletion in the gate. 

Fig. [| Drain current versus gate voltage for Vd = 1 V. Quantum mechanical treatment 
of the polysilicon gate (Ql q-poly) results in much higher current. 

Fig. |: Transmission (+) and density of states (DOS) versus energy at a spatial location 
close to the source injection barrier, at V g = OV and Vd = IV. The peaks in the density of 
states represent the resonant levels in the channel. Inset: DOS at three different y-locations 
and the total transmission. The points y = -7 and nm are to the left and right of the 
location where the source injection barrier is largest (close to y = -4 nm). 

Fig. |5|: Location of the first resonant level (E r i) versus gate voltage and the classical 
source injection barrier {E^classical)). Note that E r \ decreases slower than E^classical) 
with gate voltage due to narrowing of channel potential well. 

Fig. |6|: Plot of drain current versus gate voltage from the quantum mechanical calcula- 
tions and Medici, at Vd = IV. At small gate voltages, the drain current from Mediciil are 
comparable to the 'Ql flat band' results. The drain current from 'Ql q-poly' is however 
significantly different at large gate voltages. 

Fig. [7]: Plot of drain current versus drain voltage (V^) from the quantum mechanical 
calculations and Medici, at V g = IV. Note the large difference in drive current and dld/dVd 
between Medici,i! 'Ql flat band' and 'Ql q-poly'. 

Fig. Same as Fig. [| but the anisotropic effective mass case is included. Note that the 
valley with the largest mass in the x-direction has subband energies that are about 50 meV 
smaller than the isotropic effective mass case even at V g = 0. 

Fig. H Plot of drain current versus gate voltage for the isotropic and anisotropic effective 

1-25 



mass cases, at Vd = IV. The much higher current in the anisotropic effective mass case (Q3) 
is due to the lower suband energy shown in Fig. |8]. 

Fig. [II]: Plot of gate leakage current when the device is OFF (V g = OV) as a function of 
the y-direction, from the source to drain, for L g equal to (a) 50 nm and (b) 25 nm. Note 
the significant gate leakage current in the regions where the high doping in the source and 
drain overlap the gate in (a). A shorter gate eliminates a large fraction of the gate leakage 
current as shown in (b). 

Fig. [ll]: Polysilicon gate and oxide configurations that could reduce the OFF current 
(V g = 0V) significantly without drastically reducing the drive current (V g = IV). The 
hatched marks represent the oxide. 



Fig. |T^: Plot of drain and gate currents when the device is OFF (V g = 0V) versus oxide 
thickness for L g equal to 50 and 25 nm. Inset: Drain current for the the gate lengths when 
the device is on (V g = IV). At the larger values of t ox , the gate current (I g ) is significantly 
smaller than the drain current (Id) , meaning that the drain current is determined by electron 
injected from the source to drain. At smaller values of t ox , the drain current is dominated 
by the gate leakage current as can be seen by comparing Id and I g in this figure. More 
importantly, note that the shorter gate length (L g = 25 nm) gives an order of magnitude 
smaller drain current when the device is OFF for the smaller vaues of t ox . The inset shows 
that the drive current (Vd = V g = 1 V) is however not affected much by the shorter gate 
length. 
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FIG. 13. This figure is a larger version of Fig. 1 
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